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Inleiding 

In dit rapport wordt verslag gedaan van de eerste fase van de "GNSS Processing Groningen"^ 
uitgevoerd door de Technische Universiteit Delft (TU Delft) in opdracht van Staatstoezicht op de 
Mijnen (SodM). 

Fase 1 van de werkzaamheden betreft een analyse van de GNSS tijdreeksen van Groningen zoals 
berekend door 06-GPS. Aspecten die daarbij aan de orde kunnen komen zijn: evaluatie van de door 
06-GPS gevolgde verwerkingsmethodiek, kwantificering van de hoogteprecisie van de tijdreeksen, 
effect van de stabiliteit en invloed van de gekozen referentiestations op de tijdreeksen, effecten van 
atmosferische en andere invloeden en op de hoogtetijdseries, en correlaties tussen verschillende 
stations. 

Voor de uitvoering van de werkzaamheden is gebruik gemaakt van een GPS dataset berekend door 
het bedrijf 06-GPS in opdracht van de Nederlandse Aardolie Maatschappij (NAM). Deze dataset, met 
tijdreeksen in Latitude, Longitude en Hoogte voor 13 stations in Groningen t/m 28-02-2015, is op 25 
Maart 2015 door SodM aan de TU Delft ter beschikking gesteld. 

Voor de uitvoering van de werkzaamheden heeft de TU Delft ook toegang gekregen tot experts van 
06-GPS. Daartoe heeft op woensdag 8 April 2015 (9:00 - 11:00) overleg plaatsgevonden in Delft 
tussen Jean-Paul Henry (06-GPS) en Hans van der Marei (TU Delft). Naar aanleiding van dit overleg is 
door Frank Dentz van 06-GPS aanvullende informatie ter beschikking gesteld over de 
beschikbaarheid van referentiestations en twee performance metrics van de software. De resultaten 
hiervan zijn verwerkt in dit rapport. 

De GPS data van 06-GPS is geanalyseerd m.b.v. een aantal Matlab™ scripts voor de analyse van GNSS 
tijdsreeksen die door de auteur in de loop der jaren zijn ontwikkeld. De scripts zijn daarbij aangepast 
aan de specifieke eigenschappen van de 06-GPS data en voor de eisen van dit project. De resultaten 
worden besproken in dit rapport. 

Het eindresultaat van deze fase is dit rapport met de bevindingen en een gevalideerde tijdserie met 
kwaliteitsinformatie. Daarnaast zijn een groot aantal plots beschikbaar in een zip file. 


^ Werkzaamheden uitgevoerd volgens fase 1 van de offerte "GNSS Processing Groningen" van 27 februari en 17 
Maart 2015 en opdracht bevestiging d.d. 19 Maart 2015 met kenmerk 15039202 en inkoopordernummer 
139924. 
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GNSS dataset 


De GPS tijdsreeksen zijn door SodM beschikbaar gesteld in de vorm van twee Microsoft Excel files: 

"GPS data Groningen stations_deell2.xlsx" 

"GPS data Groningen stations_deelll.xlsx" 

De Excel files bevatten de volgende kolommen 

A number 

B doy 

C sessid (a-x) 

D number of days since 2013-01-01 00:00 (real) 

E date 

F,G,H Latitude in deg min sec 
LJ,K Longitude in deg min sec 

L Height in m 

met tab bladen voor de stations 0647, dzyl, eems, froo, over, sted (deelll) en tenp, tjuc, usqu, veen, 
zand, zdvn, zeer (deell2). 

De Matlab™ functie sodmimportdata .m (calls sodmimport stat ion .m) zet de Excel data om 
in een Matlab mat-files voor ieder station: zeer.mat zdvn.mat zand.mat veen.mat usqu.mat tjuc.mat 
tenp.mat sted.mat over.mat froo.mat eems.mat dzyl.mat 0647.mat. 


Iedere mat-file bevat een Matlab structure met de volgende velden 


'froo' 

[53.1877 6.7816 


station: 
plhO : 
epoch: 
doy: 
hour: 
plh: 
neu: 


[8736x1 

[8736x1 

[8736x1 

[8736x3 

[8736x3 


doublé] 

doublé] 

doublé] 

doublé] 

doublé] 


47.0817] <- degrees and meters 

<- Matlab datenumber 


<- degrees and meters 
<- North,East,Up ([m] wrt plhO) 


Het veld plhO bevat de gemiddelde latitude, longitude en hoogte, het veld neu bevat verschillen 
ten opzichte van deze positie in North, East en Up richting (in meters), en is berekend uit de gegeven 
latitude, longitude en hoogte gegevens. Een schatting is ieder uur beschikbaar. 

De neu tijdreeksen zijn voor een eerste inspectie geplot m.b.v. de functie sodmplotseries .m 

h=sodmplotseries('neu'^10) 

set(gca(h(l))^ 'YTick'^ [0:10:150]) 

set(gca(h(2)),'YTick',[0:10:150]) 

set(gca(h(3))^ 'YTick'^ [0:10:150]) 

print(h (1)^ '-dpng ', '-r300'^ 'SodM-North.png') 

print(h(2)^ '-dpng \ '-r300'^ 'SodM-East.png') 

print(h(3)^ '-dpng \ '-r300'^ 'SodM-Up.png') 

Het resultaat is gegeven in de onderstaande figuren. 
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Opvallend is dat de Up component veel gladder is dan zowel de North als East component, en dat 
van de twee laatsten, de North component gladder is dan de East component. Dit is consistent met 
een zogenaamde ambiguityfloat oplossing, waarbij de fase ambiguïteiten niet op hun integer 
waarden worden gefixed. Waarschijnlijker is dat de SSRPOST software van Geo++, die door 06-GPS 
wordt gebruikt, wel ambiguity fixing doet, maar dat andere parameters in de state vector 
verantwoordelijk zijn voor dit gedrag. 

Verder valt op dat sommige stations meer hoog frequente fluctuaties voortonen (ruis) dan andere 
stations, en dat vaak een duidelijk annual signaal zichtbaar is. Dit is zeker niet ongebruikelijk voor 
GPS, en we komen hier later uitgebreid op terug. Al dit soort signalen zullen we apart gaan schatten. 

Een beetje verscholen, maar toch duidelijk zichtbaar, is dat sommige stations, bv zeer, korte periodes 
hebben waarbij de data schijnbaar niet veranderd. Navraag bij 06-GPS leerde dat voor deze perioden 
waarschijnlijk geen GPS data beschikbaar is en waardoor de positie state niet wordt geupdate. Deze 
data staat weliswaar in de data files, maar het is absoluut noodzakelijk deze te verwijderen in 
verband met de verdere analyses. 

Voor het verwijderen van deze data is een Matlab functie sodmnodata. m geschreven. Deze moet 
worden aangeroepen iedere keer dat een van de mat-files wordt ingelezen (de resultaten worden in 
een andere mat file weggeschreven bij de verdere analyse van de data). Sodmnodata. m detecteert 
perioden dat de neu componenten schijnbaar niet veranderen en maakt een vector aan met de te 
verwijderen data; het is een vrij ingewikkeld script omdat er vaak valse alarm is. Uiteindelijk worden 
alleen perioden verwijderd die langer dan 18 uur duren. In alle volgende grafieken is deze data 
verwijderd. 
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In de dataset ontbreekt informatie over de kwaliteit van de individuele metingen. Het is gebruikelijk 
om voor ieder data punt een standaard afwijking afkomstig uit eerdere GPS analyse te geven. De 
waarden uit SSRPOST zijn echter veel te optimistisch en 06-GPS heeft besloten deze niet te 
verstrekken om geen valse verwachtingen te wekken. Bovendien zijn deze waarden, volgens 06-GPS, 
zeer uniform. 

Dit laatste komt ook door de gevolgde verwerking procedure door 06-GPS. De tijdseries worden 
maandelijks geupdate, door een periode van twee weken voorafgaand aan de maand (inslinger 
periode), en de maand zelf te werken. De resultaten van de laatste maand worden aan de resultaten 
van de voorafgaande berekeningen geplakt. 06-GPS doet controles op de continuïteit, en ook in de 
TU Delft analyse zijn geen verdachte sprongen op de maand overgangen gevonden. 


GNSS data fitting procedure 

GNSS tijdseries worden door drie verschillende effecten beïnvloed: 

1. Beweging van het station zelf (in een welbepaald referentieframe). 

2. Beweging van het monument onder invloed van de omgeving, aardgetijden, atmospheric en 
ocean loading, grondwater effecten, etc.. Deze bewegingen worden deels gecorrigeerd in 
de GNSS data analyse, waaronder met name de getijde en loading effecten door hetzij het 
toepassen van expliciete correcties, of anderzijds door het vasthouden van 
referentiestations (netwerk processing). De impliciete aanname in de netwerk processing is 
dat de effecten voor alle stations gelijk zijn; dit is echter slechts deels waar. Bv voor 
Groningen, ocean loading is niet hetzelfde voor alle stations en het is niet duidelijk of de 
SRRPOST software voor dit effect corrigeert. Het gevolg is vaak dat periodieke effecten in de 
GNSS tijdseries zichtbaar zijn, als direct effect, of door aliasing met een gekozen data 
interval. 

3. Schijnbare effecten t.g.v. bijvoorbeeld niet gemodelleerde antenne effecten (calibraties), 
site multipath, niet gemodelleerde atmosfeer, etc. Correcties in de GNSS processing zijn 
vaak elevatie (en azimuth) afhankelijk en deze kunnen doorwerken in de tijdseries onder 
invloed van de zich herhalende satelliet constellatie. 

De twee laatst genoemde effecten werken verstorend voor de analyse van bodembeweging in 
Groningen. Ze zullen proberen deze uit te data te schatten en vervolgens de tijdseries hiervoor 
corrigeren, voordat we een analyse van de bodembeweging kunnen uitvoeren. 

Stations model 

Iedere component AN, AE, AU (North, East, Up) van de positie tijdseries wordt beschreven door het 
volgende model 


A = s(t) -I- AMrnLdiP - Po) + ^Tempi(J “ W + T,(asiSin2nfit -I- UciCosInfit ) -I- e 

met s(t) de trend, een coëfficiënt voor atmospheric loading en A^g^pj een coëfficiënt voor 

het temperatuureffect, Ogi en a^i coëfficiënten voor de harmonische termen en e de residuele 
noise, t de tijd in decimal years en fi de frequentie in cycles/year. Het model heeft voorts nog 
mogelijkheden voor het schatten van jumps, maar deze is in dit project niet nodig gebleken. 


5 


Normaliter wordt een jump gemodelleerd wanneer er veranderingen in de (antenne) hardware 
plaatsvinden, of anderszins bij sprongen in de data. 

Het trend model s(t) kan een simpele lineaire trend zijn, maar ook hogere orde polynomen en 
(cubic) splines zijn mogelijk. Omdat de meeste tijdseries slecht een jaar lang zijn wordt begonnen 
met een lineair model s(t) = x(tO) + v * (t — tO). Alleen voor tenp wordt in een later stadium een 
cubic spline model gebruikt. 

In het standaard model wordt een correctie geschat voor atmospheric loading en stations 
deformaties onder invloed van temperatuur. Daarbij wordt gebruik gemaakt van de luchtdruk P en 
temperatuur T van naburige meteo stations, in dit geval het KNMI weerstation in Eelde. Vanwege de 
beperkte omvang van het gebied zal het atmospheric loading effect zal heel erg klein zijn, maar 
wordt desalnietemin meegeschat. De geschatte coëfficiënten zijn ook daadwerkelijk klein en we 
hadden de schatting van dit effect ook achterwege mogen laten, maar het meenemen van deze term 
schaad ook niet. 

De harmonische termen die worden geschat zijn periodes van 1 cycle/jaar en 2.08 cycle/jaar: een 
annual en semi-annual periode. De semi-annual periode is 175 dagen, de helft van de GPS draconitic 
jaar van 351 dagen en de helft van de periode van 350 dagen dat de satelliet constellatie zich 
herhaalt. De GPS draconitic jaar periode zelf wordt niet meegeschat aangezien deze te dicht bij de 
annual periode zit. Deze twee perioden zijn zeer gebruikelijk in de analyse van GPS tijdseries. Het is 
mogelijk dat ook hogere harmonischen van de GPS draconitic periode voorkomen, maar uit de 
verdere data analyse zal blijken dat het niet nodig is deze er bij te doen. 

De temperatuur vertoond ook een duidelijke jaarlijkse oscillatie. Toch is het mogelijk het 
temperatuur effect en de annual terms gezamenlijk te schatten zonder dat de individuele 
schattingen slecht worden. Dit komt doordat de temperatuur veranderingen van jaar tot jaar 
verschillen en de temperatuur veranderingen geen puur sinusvormig patroon volgen. 

Gecorrigeerde tijdserie 

De verschillende parameters worden geschat met behulp van de kleinste kwadraten methode. Nadat 
de parameters zijn geschat worden de kleinste kwadraten residuen ê berekend. De gecorrigeerde 
tijdserie, zonder periodieke of verstorende effecten, is 

A = s(t) + ê 


met s(t) de geschatte trend. 

Meteo data KNMI station Eelde (280) 

Het Matlab script gmeteo. m leest uurlijkse meteo data van het KNMI weerstation Eelde. Het script 
berekend ook het dagelijkse gemiddelde van de luchtdruk en temperatuur, en plot deze. De 
resultaten voor de temperatuur staan in de onderstaande figuren. Deze plots, samen met die voor de 
luchtdruk (niet afgedrukt), zijn ook beschikbaar in een zip file (Eelde_Yearly_Temp.eps, 
Eelde_Yearly_Temp.png, Eelde_Yearly_Pressure.eps, Eelde_Yearly_Pressure.png, Eelde_Temp.png, 
Eelde_Temp.eps. Eelde_Pressure.png, Eelde_Pressure.eps). De data wordt opgeslagen in een mat file 
zlmeteo .mat die door sodmtseriesplotexp .m wordt gelezen. 
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Voor het schatten van de temperatuur invloed wordt gebruik gemaakt van een interpolatie op basis 
van de dagelijkse gemiddelden (de uurlijkse meteo is niet geschikt ivm de grote dag/nacht variaties). 


GNSS data fitting resultaten - Eerste iteratie 

De hierboven beschreven GNSS data fitting procedure is geïmplementeerd in de Matlab functie 

tseriesanalysis .m. Deze functie wordt aangeroepen vanuit de functie 

sodmtseriesplot. m. De functie sodmtseriesplot. m voert de volgende taken uit: 

leest de mat files met stations data, 

detecteerd m.b.v. sodmnodata .m intervallen met ongeldige metingen, 
leest de meteo data uit de mat file zlmeteo .mat, 
berekend decimal year, 

voert de timeseries analysis uit door een aanroep van tseriesanalysis .m, 
schrijft de resultaten van de timeseries analysis weg naar een mat-file met de naam 
ssss_f it .mat, waarbij ssss de stations naam is, 
plot de resultaten. 

Voor ieder station worden drie plots geproduceerd: ssss_series .png met de raw data series en 
fit, ssss_components . png met de geschatte componenten en ssss_residuals . png met 
de residuals. Voorbeelden van de plots voor station tenp worden hieronder gegeven. 
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De plots voor de overige stations zijn beschikbaar in de zip file met plots. 
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De cyaan achtergrond is het interval [-lai ZcrJ , met (Tj de standaard afwijking van de individuele 
waarnemingen in de tijdreeks. Gebruikelijk is deze standaard afwijking afkomstig is uit de eerdere 
GPS analyse, echter aangezien deze ontbreekt in de dataset is een waarde van 1 mm aangenomen. 
De waarden uit SSRPOST zijn te optimistisch en 06-GPS heeft besloten deze niet te verstrekken om 
geen valse verwachtingen te wekken. Wel zijn deze waarden, volgens 06-GPS, zeer uniform. Dit is 
een belangrijk gegeven aangezien de absolute grootte van niet van belang is voor de schattingen 
(alleen voor de kwaliteitsbeschrijving). Ook wordt tijdens de processing een schalingsfactor (OMT) 
geschat. 

Om de interpretatie te vergemakkelijken zijn de resultaten ook voor alle stations tegelijk geplot, 
vergelijkbaar met de figuur op pagina 3 en 4. Hiervoor is gebruik gemaakt van de functie 

sodmplotseries .m en het script sodmplotsum .m. 

De resultaten voor de verschillende componenten worden hieronder gegeven: 

de originele GNSS data met ongeldige data verwijderd (raw) 
de model fit (fit) 

de som van de harmonische termen en de temperatuur invloed (harmonie + tempi) 
de residuen na de fit 

Deze plots en een aantal anderen zitten ook in de zip file met plots onder de naam all_* . png. 


Groningen (raw) 
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De plots laten belangrijke verschillen zien tussen de individuele stations, maar ook een aantal 
overeenkomsten. 

Voor de validatie van de model fit is een nadere inspectie van de residuen belangrijk. Een aantal 
zaken valt daarbij op 

de belangrijkste periodieke effecten zijn verdwenen, maar de residuen zijn duidelijk 
gecorreleerd in de tijd en bevatten mogelijk nog hogere harmonischen of andere effecten 
alle residual series lijken op het eerste gezicht gemiddeld om de nul waarde te schommelen 
(zero mean), met uitzondering van de hoogte van tenp dat alle schijn heeft dat het linaire 
model niet goed fit, en met uitzondering van een korte periode aan het begin in 2013 waarbij 
het lijkt of een sprong zit, 

de residuen in de East component het grootst zijn, 
dat er correlatie is tussen stations onderling. 


Later zullen we een meer gedetaileerde analyse van de residuen uitvoeren. 


De geschatte parameters worden in de onderstaande tabel gegeven: 




Vel 

AtmLd 

Templ 



mm/y 

mm/kPa 

mm/daK 

0647 

Lat 

-1.76 

-0.00 

-0.09 

0647 

Lon 

3.95 

0.03 

-1.82 

0647 

Rad 

-1.40 

0.04 

-0.58 

dzyl 

Lat 

-0.10 

0.15 

0.18 

dzyl 

Lon 

1.33 

-0.07 

-0.08 

dzyl 

Rad 

-2.71 

-0.03 

-0.39 


365d 

175d 

StdF 

StdR 

OMT 

mm 

mm 

mm 

mm 


2.16 

0.19 

1.00 

0.65 

0.42 

2.44 

0.39 

1.00 

1.08 

1.16 

3.03 

0.37 

1.00 

0.82 

0.66 

0.76 

0.47 

1.00 

0.63 

0.39 

0.63 

0.34 

1.00 

0.85 

0.73 

1.96 

0.12 

1.00 

0.66 

0.44 
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eems 

Lat 

-1.72 

-0.02 

-2.15 

1.62 

0.12 

1.00 

0.67 

0.45 

eems 

Lon 

-2.22 

-0.04 

0.85 

0.34 

0.09 

1.00 

0.75 

0.56 

eems 

Rad 

-5.33 

-0.00 

-0.37 

0.74 

0.64 

1.00 

0.55 

0.30 

f roo 

Lat 

0.76 

-0.04 

-0.23 

0.51 

0.40 

1.00 

0.47 

0.22 

f roo 

Lon 

0.88 

-0.04 

2.37 

2.35 

0.74 

1.00 

0.97 

0.94 

f roo 

Rad 

-5.76 

-0.14 

-0.06 

1.16 

0.69 

1.00 

0.64 

0.41 

over 

Lat 

-2.07 

0.06 

-2.35 

1.80 

0.28 

1.00 

0.72 

0.52 

over 

Lon 

-2.35 

-0.00 

-0.99 

1.33 

0.07 

1.00 

0.80 

0.65 

over 

Rad 

-2.66 

0.01 

-0.18 

1.52 

0.10 

1.00 

0.63 

0.39 

sted 

Lat 

-0.32 

0.03 

0.00 

0.50 

0.19 

1.00 

0.52 

0.27 

sted 
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0.10 

0.01 

-1.01 

1.06 

0.21 

1.00 

0.77 

0.59 

sted 

Rad 

-4.76 

-0.01 

-0.43 

0.42 

0.36 

1.00 

0.53 

0.28 

tenp 

Lat 

-1.32 

0.12 

1.98 

1.25 

0.37 

1.00 

0.91 

0.83 

tenp 

Lon 

0.58 

0.18 

1.41 

1.15 

0.31 

1.00 

1.24 

1.55 

tenp 

Rad 

-7.88 

0.02 

-0.94 

1.94 

0.34 

1.00 

0.92 

0.85 

t j uc 

Lat 

-0.19 

0.02 

-0.32 

0.44 

0.14 

1.00 

0.50 

0.25 

t j uc 

Lon 

3.05 

0.01 

-0.06 

1.01 

0.32 

1.00 

0.67 

0.46 

t j uc 

Rad 

-3.82 

0.02 

-0.26 

1.24 

0.25 

1.00 

0.71 

0.50 

usqu 

Lat 

-1.76 

0.06 

0.17 

0.78 

0.36 

1.00 

0.55 

0.31 

usqu 

Lon 

-0.24 

-0.02 

-0.72 

0.72 

0.40 

1.00 

0.74 

0.55 

usqu 

Rad 

-2.65 

-0.01 

-0.56 

0.15 

0.48 

1.00 

0.67 

0.45 

veen 

Lat 

3.33 

0.01 

0.22 

0.65 

0.22 

1.00 

0.75 

0.56 

veen 

Lon 

7.14 

0.09 

0.97 

0.48 

0.27 

1.00 

0.93 

0.86 

veen 

Rad 

-5.96 

-0.20 

-0.36 

1.35 

0.36 

1.00 

0.84 

0.70 

zand 

Lat 

-1.26 

-0.06 

-2.01 

1.92 

0.05 

1.00 

0.67 

0.44 

zand 

Lon 

-4.41 

-0.08 

-2.42 

1.91 

0.32 

1.00 

0.96 

0.92 

zand 

Rad 

-4.18 

0.07 

-0.13 

0.44 

0.26 

1.00 

0.70 

0.49 

zdvn 

Lat 

0.04 

-0.05 

0.36 

0.36 

0.30 

1.00 

0.47 

0.22 

zdvn 

Lon 

-0.83 

0.09 

3.52 

1.89 

1.20 

1.00 

1.19 

1.41 

zdvn 

Rad 

-4.68 

-0.08 

0.16 

1.14 

0.57 

1.00 

0.69 

0.47 

zeer 

Lat 

-1.53 

-0.05 

-0.32 

1.13 

0.71 

1.00 

0.66 

0.44 

zeer 

Lon 

2.20 

0.03 

-0.22 

0.83 

0.47 

1.00 

0.77 

0.59 

zeer 

Rad 

-5.28 

0.11 

-0.52 

0.82 

0.31 

1.00 

0.56 

0.31 


Voor de harmonische termen is alleen de amplitude afgedrukt. De precisie van de geschatte 
parameters beter dan 0.1 mm. 

StdF is het gemiddelde van de Standard afwijking die gegeven is voor de waarnemingen (in ons geval 
de gekozen waarde van 1 mm). StdR is de emperische standaard afwijking van de residuen en OMT 
zijn de resultaten van de Overall Model Toets 


OMT = 



Aangezien de standaard afwijking voor alle waarnemingen gelijk is (1 mm) is deyfOMT vrijwel 
gelijk aan de empirische standaard afwijking van de residuen (het enige verschil is de noemer in de 
deling). 


Er zijn grote verschillen zichtbaar tussen de stations voor de harmonische componenten en de 
geschatte temperatuur invloed. Echter er zijn ook overeenkomsten tussen sommige stations. Daarom 
zijn de belangrijkste parameters ook geplot in kaart vorm m.b.v. de functie socimplotmap . m . De 
oorzaken en verbanden tussen de harmonischen en de temperatuur invloed zijn niet verder 
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onderzocht. Oorzaken zijn mogelijke verschillen in constructie van de stations (andere gebouwen), 
maar er kan ook gekeken worden naar residuele effecten van ocean loading. De atmospheric loading 
term is zoals verwacht verwaarloosbaar klein voor alle stations en is niet apart geplot. 


Groningen GPS (Annual Term) Groningen GPS (Temp.Influence) 



Groningen GPS (Velocity) Groningen GPS (Emperical co-variance) 



De plot rechtsonder toont de empirische co-variantie zoals berekend uit de residuen. Ook hier valt op 
dat de empirische covariantie van tenp groter is dan de overige stations wat op een slechte model fit 
duidt. 
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GNSS Spectra 

Met behulp van de Matlab functie sodmspectra. m wordt de Lomb-Scargle periodogram 
berekend voor verschillende grootheden. Dit kan voor individuele stations, maar het is tevens 
mogelijk de periodogrammen van alle stations te stapelen (stacken) voor een betere resolutie. 

In de onderstaande figuren wordt het stacked periodogram gegeven van het detrended signaal 
(alleen lineaire trend) en het stacked periodogram van de residuen (dus na verwijdering van een 
lineaire trend, temperatuur invloed, annual en semi-annual termen). 




De North, East en Up componenten zijn voor de leesbaarheid ten opzichte van elkaar verschoven: de 
gestippelde schuine lijnen geven het referentie niveau van 0.1 f~^ [mrn}/cpy] (flicker noise). De 
vertikale gestippelde lijnen vertegenwoordigen harmonische perioden typisch voor GPS: de annual 
periode van 1 cpy en semi-annual van 0.5 cpy, de GPS draconitic periode van 1.04 cpy en de hogere 
harmonischen daarvan, 2.08, 3.12, 4.16, 5.20 en 6.24 cpy, en een periode van 14.75 dagen. 
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Behalve de verwijderde 1 cpy en 2.08 cpy perioden zijn er geen overtuigende andere perioden die op 
dit moment geschat moeten worden. In de Up series zit nog wel een duidelijke piek bij 7.28 cpy (de 
T harmonische van de draconitic period) en in de North component bij 13.22 cpy (27.63 days) en 
3.13 cpy in het stacked periodogram, maar deze zijn minder goed waarneembaar in de individuele 
periodogrammen. Het is mogelijk dat deze bij residual stacking (zie volgende sectie) verdwijnen. 

De signalen vertonen een typische power law gedrag overeenkomstig flicker noise. De up 
component vertoont zelfs lichte trekjes van een random walk (ƒ “^), maar is toch nog steeds 
overwegend flicker noise. 


Variantie component schatting 

De periodogrammen stellen ons in staat het kansmodel voor de GPS waarnemingen te bepalen. In de 
plaats van ongecorreleerde white noise met en standaardafwijking van 1 mm, hetgeen resulteert in 
een diagonale covariantie matrix Qy voor de waarnemingen, is het realistischer met een volle 
covariantie matrix gegenereerd op basis van de aanname van flicker noise te werken. Zonder deze 
zijn alle schattingen van de varianties te optimistisch. De schalings factoren zijn uit de OMT test te 
bepalen. Nog een andere oplossing is variantie component schatting toe te passen; m.b.v. variantie 
component schatting is het mogelijk het aandeel white noise, flicker noise en brownian (random 
walk) noise te bepalen, of de macht in de power law. 

Een verbeterd kansmodel draagt vervolgens weer bij in verbeterde schattingen van de snelheden, 
temperatuurs effect, harmonische termen, etc., en realistische schattingen van de standaard 
afwijking van deze parameters. Deze stap is, vanwege de beperkte tijd, in dit rapport nog niet 
uitgevoerd , maar wordt wel aanbevolen om dit te gaan doen om zodoende realistische schattingen 
te krijgen van de standaard afwijking van verschillende parameters, waaronder de snelheden. 


Common mode - Residual stack 

Om mogelijke common mode signalen te onderzoeken zijn de residuals gestapeld (stacked). De 
Matlab functie sodmresidualstack .m berekent het gemiddelde en de standaard afwijking van 
de residuen over alle stations. Het gemiddelde en de standaardafwijking kan naar keuze berekende 
worden over het inteval van een uur of een dag. 

De onderstaande plot laat het uurlijkse gemiddelde, het dagelijkse gemiddelde en uurlijkse standaard 
afwijking zien. In een tweede plot zijn het dagelijkse gemiddelde en standaard afwijking geplot. Het 
dagelijkse gemiddelde is dezelfde orde grootte, en soms groter, als de dagelijkse standaardafwijking. 
De standaard afwijking van het gemiddelde zelf is lager, met een factor 1.7 voor de begin periode 
met 3 stations, en een factor 3.6 voor de periode na Feb 2014 met 13 stations. De voorzichtige 
conclusie is dat het gemiddelde niet zo maar ruis is, maar een common mode signaal 
vertegenwoordigd. Het heeft daarom zin om in een volgende iteratie van de fitting procedure te 
corrigeren voor het common mode signaal uit residual stacking. 
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A U [mm] A E [mm] A N [mm] A U [mm] A E [mm] A N [mm] 


Groningen (Stacked Residuals) 





Mean (hourly) Mean (daily) Std (hourly) 


Groningen (Stacked Residuals) 





Mean (daily) Std (daily) 














De residuals zijn als test ook per dag van de maand gestacked. Dit heeft te maken met het gegeven 
dat 06-GPS de processing maandelijks opstart, te beginnen met een inslinger periode van twee 
weken, en gevolgd door een maand data die vervolgens aan de eerdere tijdserie wordt geplakt. De 
resultaten zijn hieronder weergegeven en laten zien dat deze procedure werkt en geen nadelige 
effecten in de tijdserie geeft. 


Groningen (Stacked Residuals) - Monthiy 





Mean Std 


De Stacked residuals zijn niet alleen te gebruiken als correcties voor een volgende iteratie, maar het 
heeft ook zin om te controleren of deze niet correleren met andere (performance) parameters uit de 
SSRPOST processing en de beschikbaarheid van data van de gebruikte referentiestations. 

Om dit mogelijk te maken heeft 06-GPS een aantal performance parameters beschikbaar gesteld die 
in de volgende secties worden geanalyseerd. 


Irregularity parameters (IR, IPi en IPO) uit de SSRPOST processing 

SSRPOST bewaard ionosfeer/troposfeer log-bestanden met de irregularity parameters IR, IPi en IPO 
De irregularity parameters worden iedere minuut opnieuw berekend. De waarden zijn beschikbaar 
van zowel de hoofdserver (slie) als van de back-up server (deve) beschikbaar. Er kunnen er soms 
valse uitschieters in de data zitten. Daarom zijn voor ieder uur en iedere dag de 95% percentilen 
berekend. Deze dagelijkse waarden zijn geplot in de onderstaande figuur. De uurlijkse waarden zijn 
beschikbaar in de zip file met plots. 
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IR, IPi and IPO daily 95% values deve 



De volgende beschrijving van de parameters is gecopieerd uit de handleiding van SSRPOST: 

"The value IR is computed from the Standard deviation of the stochastic 3D Gauss-Markov processfor one satellite. The IR 
number shown here is the maximum Standard deviation currently present in the system from all satellites. The IR irregularity 
is like the 195 also not suited os on indicator of the service quolityfrom o rover user's point of view, but con give the RTK 
network provider useful Information on the maximum ionospheric process Standard deviation. 

The two IR Irregularity values ore o meosure of the interpolotion quolity of the state Information, i.e. how well the network 
con predict the distance dependent errors in between the reference stations. They ore disployed (in meter) separoteiy for the 
ionospheric (IP-I) and the geometrie (IP-0) port of the error components. They ore computed os the Standard deviation of o 
second order FKP representotion of the current state spoce using the Standard deviation computed from all stations and 
satellites with o distance dependent weighting. The IP reodings correspond to the remoining third or higher order effects. 

The numbers are derived from the non-describoble residuol errors with the ossumption of o quodrotic FKP model (which 
must not necessorily exoctly correspond to the models used in the octuol RTK network). 

A RTK network user equipment hos to expect and must account for these IP residuol errors in the field. The ionospheric port 
IP-I con normolly be eliminoted by duol frequency receivers, ifthe rover could solve lts ombiguities. The geometrie port IP-0 
is moiniy influenced by tropospheric irregulorities (since satellite orbits today normolly don't show significant irreguiorities 
and have thus in proctice no influence on IP-0). Unfortunotely no Standard dato formot is aval la be today to teil the rover 
these IP residuols. It is however possible through the distance of the Pseudo Reference Station PRS (see parameter o in the 
RTCM_OUT option -FKP,V,o) to give to the rover on indirect hint, whot magnitude of residuols con be possible. 

The reliable computation IP Irregularity values requires o large redundancy in the network, while the IR Irregularity values 
alreody ollow on interpretation even for low redundoncies. In case of low redundancy, afirst order FKP representotion (i.e. 
no quodrotic FKP option -q) should be used for the computation of the irregulorities. Ifthere is no redundancy at all (e.g. too 
few stations) then no IP Irregularity values ore computed and the fields IP-I and IP-0 remoin empty. 


For o detail led theoreticol discussion of the 195 parameters and the Irregularity parameters refer to the Geo++® White 
Paper: Gerhord Wübbeno, Martin Schmitz, Andreos Bogge: GNSMART Irregularity Reodings for Distance Dependent Errors 
Download as PDF - size 203 kB." 
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Op het eerste gezicht lijkt een lichte correlatie te bestaan tussen de Standard afwijking van de 
stacked residuals en de IPO (troposfeer) parameter. Dit soort aspecten zou echter nog verder 
onderzocht kunnen worden. 


Missing data for the reference stations 

De ontbrekende data van de referentiestations is door 06-GPS samengevat in de excel file 
missing_epochs_NAM. xlsx . Deze informatie is echter nog niet geanalyseerd. 


GNSS data fitting resultaten - Tweede iteratie 

De tweede tweede iteratie van de data fitting procedure verschilt van de eerste in de volgende 
aspecten: 

de data is corrigeerd voor de common mode signal uit de stacked residuals van de eerste 
iteratie, 

voor Ten Post (tenp) is een cubic spline model bestaande uit twee piecewise polynomial 
intervallen en continuïteit in de 1® afgeleide tussen de piecewise polynomials. 

Het kansmodel is nog niet aangepast. Geadviseerd wordt om i.p.v. een white noise model een flicker 
noise model te gebruiken en variantie component schatting toe te passen. Dit zou een strenge 
toetsing van de resultaten mogelijk maken en resulteren in realistische schattingen voor de 
standaard afwijking van de parameters. 

De tweede iteratie is geïmplementeerd in de functie sodmtseriesplotexp .m en script 
sodmtseriesiter .m. Voor ieder station worden drie plots geproduceerd: 

ssss_series2 .png met de raw data series en fit, ssss_components2 .png met de geschatte 
componenten en ssss_residuals2 .png met de residuals. 

Daarnaast worden een aantal overzichtplots gemaakt van de gecorrigeerde tijdseries en residuen. De 
gecorrigeerde tijdseries zijn hieronder geplot. De annual term, semi-annual term, temperatuur 
invloed, atmospheric loading en common mode signalen zijn verwijderd. De gecorrigeerde tijdserie is 
berekend uit de geschatte trend s(t) en residuen ê 

A = ■s(t) + 6 , 

waarbij de trend een linear model is voor alle stations met uitzondering van tenp. 

De overige plots zijn beschikbaar in de zip file met plots. 
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A E [mm] A U [mm] 


Groningen (trend + residual) 
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Groningen (trend + residual) 
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De geschatte parameters zijn in de onderstaande tabel gegeven: 




Vel 

AtmLd 

Templ 

365d 

175d 

StdF 

StdR 

OMT 



mm/y 

mm/kPa 

mm/daK 

mm 

mm 

mm 

mm 


0647 

Lat 

-1.80 

0.02 

-0.14 

2.15 

0.13 

1.00 

0.48 

0.23 

0647 

Lon 

3.88 

0.02 

-1.85 

2.50 

0.34 

1.00 

0.89 

0.79 

0647 

Rad 

-1.36 

0.04 

-0.66 

3.10 

0.23 

1.00 

0.66 

0.43 

dzyl 

Lat 

-0.47 

0.13 

0.21 

0.87 

0.41 

1.00 

0.46 

0.21 

dzyl 

Lon 

1.04 

-0.04 

0.01 

0.49 

0.30 

1.00 

0.69 

0.47 

dzyl 

Rad 

-3.40 

-0.03 

-0.41 

1.84 

0.09 

1.00 

0.47 

0.22 

eems 

Lat 

-1.77 

-0.03 

-2.12 

1.66 

0.16 

1.00 

0.55 

0.30 

eems 

Lon 

-1.76 

-0.03 

0.90 

0.24 

0.04 

1.00 

0.58 

0.34 

eems 

Rad 

-5.56 

-0.00 

-0.42 

0.72 

0.61 

1.00 

0.42 

0.18 

f roo 

Lat 

0.62 

-0.05 

-0.19 

0.58 

0.38 

1.00 

0.39 

0.16 

f roo 

Lon 

1.18 

-0.04 

2.42 

2.37 

0.69 

1.00 

0.78 

0.60 

f roo 

Rad 

-6.10 

-0.14 

-0.10 

1.16 

0.66 

1.00 

0.47 

0.22 

over 

Lat 

-2.14 

0.05 

-2.31 

1.84 

0.27 

1.00 

0.55 

0.30 

over 

Lon 

-1.96 

0.01 

-0.94 

1.29 

0.15 

1.00 

0.63 

0.40 

over 

Rad 

-2.93 

0.01 

-0.22 

1.50 

0.14 

1.00 

0.38 

0.15 

sted 

Lat 

-0.84 

0.02 

0.08 

0.68 

0.12 

1.00 

0.33 

0.11 

sted 

Lon 

-0.08 

0.03 

-0.90 

0.97 

0.18 

1.00 

0.55 

0.31 

sted 

Rad 

-5.12 

-0.01 

-0.40 

0.31 

0.35 

1.00 

0.42 

0.18 

tenp 

Lat 

-1.89 

0.11 

2.01 

1.20 

0.32 

1.00 

0.73 

0.53 

tenp 

Lon 

-0.44 

0.12 

1.59 

1.20 

0.33 

1.00 

0.94 

0.87 

tenp 

Rad 

-8.65 

-0.02 

-0.40 

1.62 

0.41 

1.00 

0.51 

0.26 

tjuc 

Lat 

-0.23 

0.01 

-0.28 

0.49 

0.17 

1.00 

0.32 

0.10 

tjuc 

Lon 

3.48 

0.02 

-0.02 

1.09 

0.29 

1.00 

0.51 

0.26 

tjuc 

Rad 

-4.07 

0.01 

-0.31 

1.27 

0.26 

1.00 

0.49 

0.24 
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usqu 

Lat 

-2.36 

0.04 

0.20 

0.99 

0.28 

1.00 

0.45 

0.20 

usqu 

Lon 

-0.36 

0.01 

-0.65 

0.65 

0.39 

1.00 

0.64 

0.40 

usqu 

Rad 

-2.99 

-0.00 

-0.55 

0.25 

0.45 

1.00 

0.63 

0.39 

veen 

Lat 

3.29 

0.03 

0.18 

0.63 

0.15 

1.00 

0.56 

0.32 

veen 

Lon 

7.07 

0.08 

0.94 

0.42 

0.22 

1.00 

0.73 

0.54 

veen 

Rad 

-5.92 

-0.20 

-0.45 

1.42 

0.46 

1.00 

0.70 

0.49 

zand 

Lat 

-1.30 

-0.07 

-1.98 

1.95 

0.08 

1.00 

0.54 

0.29 

zand 

Lon 

-3.97 

-0.08 

-2.38 

1.86 

0.41 

1.00 

0.87 

0.76 

zand 

Rad 

-4.43 

0.07 

-0.18 

0.44 

0.24 

1.00 

0.54 

0.29 

zdvn 

Lat 

-0.05 

-0.05 

0.40 

0.42 

0.27 

1.00 

0.38 

0.15 

zdvn 

Lon 

-0.47 

0.10 

3.57 

1.90 

1.13 

1.00 

1.03 

1.06 

zdvn 

Rad 

-4.99 

-0.08 

0.11 

1.17 

0.53 

1.00 

0.52 

0.28 

zeer 

Lat 

-1.93 

-0.06 

-0.31 

1.21 

0.77 

1.00 

0.49 

0.24 

zeer 

Lon 

2.04 

0.06 

-0.19 

0.81 

0.44 

1.00 

0.55 

0.30 

zeer 

Rad 

-5.60 

0.11 

-0.54 

0.82 

0.29 

1.00 

0.43 

0.19 


De kwaliteit van de fit is aanzienlijk verbeterd: de emperische standaardafwijking StdE en overall 
model test OMT zijn aanzienlijk lager dan in de eerste iteratie. De snelheid voor tenp in de tabel en 
plots is de snelheid in het midden van het data interval (snelheid voor tenp is niet constant). 

De belangrijkste parameters en de covariantie zijn evenals in de eerste iteratie ook als map 
weergegeven: 


Groningen GPS (Velocity) Groningen GPS (Emperical co-variance) 
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Groningen GPS (Annual Term) Groningen GPS (Temp.Influence) 



De stacked Lomb-Scargle periodogram van de residuen en de residual stack voor de tweede iteratie 
zijn hieronder weergegeven. 
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Groningen (Stacked Residuals) 



Mean (hourly) Mean (daily) Std (hourly) 


Ten Post (TENP) 

De onderstaande plot laat de uiteindelijke resultaten voor Ten Post (tenp) zien. De getoonde 
tijdseries is de gecorrigeerde tijdserie uit de tweede iteratie, zonder periodieke of verstorende 
effecten, en is berekend uit de cubic spline fit s(t) en residuen ê 

A = s(t) + ê 

De cubic spline bestaat twee piecewise polynomial intervallen en continuïteit in de Ie afgeleide 
tussen de piecewise polynomials. De spline fit en residuen zijn bepaald uit de met de residual stack 
gecorrigeerde data. 
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Corrected Series-Fit 


De daling in Ten Post laat een afnemende trend zien. 


De verschillen tussen de eerste en tweede iteratie, en tussen een lineair en spline model zijn 
hieronder weergegeven. 
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TENP 



Spline VS linear model-First vs second iteration 


Wat direct opvalt is dat het onderliggende trend model relatief weinig van invloed is aangezien deze 
wordt gecorrigeerd met de residuen. 


Conclusie 

Het doel van dit onderzoek was de GNSS tijdreeksen van Groningen zoals berekend door 06-GPS 
nader te analyseren. Hiervoor is een analyse methode opgezet die de tijdreeksen ontleed in 
verschillende signalen: een lineaire of spline bewegingsmodel, annual en een semi-annual signaal, 
temperatuur invloed, common mode signaal, en een residueel signaal. 

Door de annual, semi-annual, temperatuur invloed en common mode signalen te verwijderen uit de 
tijdreeksen ontstaat een duidelijker beeld van de eigenlijke bodembeweging: het doel van de 
metingen. De combinatie van een geschatte trend en geschatte residuen, na correctie voor een 
common mode signaal, geeft een hoge resolutie bodembewegingsignaal. De resultaten laten ook 
zien dat dit signaal relatief ongevoelig is voor het onderliggende trend model. 

De GNSS metingen op het station in Ten Post laten zien dat de daling gelijdelijk afneemt en goed kan 
worden beschreven met een spline model. 

Gezien de beperkt beschikbare tijd voor dit onderzoek zijn niet alle vragen in evenveel detail 
onderzocht. 


De belangrijkste vraag die nog open staat is het kansmodel van de GNSS metingen. Een eerste aanzet 
is gegeven in dit onderzoek en het is duidelijk dat de waarnemingen overwegend beïnvloed worden 
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doorflicker noise. Echter, in de berekeningen is nog uitgegaan van een white noise model. Dit is niet 
alleen van invloed op de geschatte parameters, maar heeft ook tot gevolg dat a) nog geen 
realistische beschrijving gegeven kan worden van de standaard afwijking van de geschatte 
parameters (bv de standaard afwijking van de geschatte snelheden of snelheidsveranderingen), en b) 
verschillende model aannames niet streng getoetst kunnen worden. 


[einde van dit document] 
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